# Turn on the gmri font for plots
showtext::showtext_auto()

Ocean warming Patterns of New England

This report seeks to showcase how different areas in New England are warming in the larger context of the rest of the world. The two main focal areas are the Gulf of Maine and the Northeastern United States’ Continental Shelf.

Regional Timeseries Patterns

Loading Timeseries

# OISST Data
gom_ts <- oisst_access_timeseries(oisst_path = box_paths$oisst_mainstays, 
                                  region_family = "gmri focus areas", 
                                  poly_name = "apershing gulf of maine")
shelf_ts <- read_csv(get_timeseries_paths("nelme_regions")$NELME$timeseries_path)
world_ts <- read_csv(str_c(oisst_path, "/global_timeseries/global_anoms_1982to2011.csv"))


# Put them in a list 
area_list <- list(
  "Global Oceans" = world_ts,
  "Northeastern U.S. Shelf" = shelf_ts,
  "Gulf of Maine" = gom_ts)

# Make time a date
area_list <- map(area_list, ~ mutate(.x, time = as.Date(time)))


# Get Anomalies and heatwave info
area_hw <- map(area_list, pull_heatwave_events)

Summarizing Data

# 3. Get monthly averages
area_means <- map_dfr(area_hw, function(x){
  x %>%
    mutate(mnth = month(time, label = T, abbr = T),
           yr = year(time)) %>% 
    group_by(yr, mnth) %>% 
    summarise(avg_temp  = round(mean(sst, na.rm = T), 2),
              avg_anom  = round(mean(sst_anom, na.rm = T), 2),
              n_hw_days = round(sum(mhw_event, na.rm = T)), 2) %>% 
    ungroup() 
  
  }, .id = "Region")

2021 Tables

The following tables are designed to highlight the relative magnitude of the local events of the Gulf of Maine or NE Shelf with the rest of the world’s oceans.

Format Data for Tables

# Pivot for Table Displays
anom_avgs <- area_means %>% 
  filter(yr == 2021) %>% 
  select(Region, mnth, avg_anom) %>% 
  pivot_wider(names_from = mnth, values_from = avg_anom)


heat_avgs <- area_means %>% 
  filter(yr == 2021) %>% 
  select(Region, mnth, n_hw_days) %>% 
  pivot_wider(names_from = mnth, values_from = n_hw_days)

Monthly Anomalies

# 1. Anomalies
anom_avgs %>% 
  gt(rowname_col = "Region") %>% 
  tab_stubhead(label = "Region") %>% 
  tab_header(
    title = md("**Average Temperature Anomalies - 2021**"), 
    subtitle = paste("Degree Celsius Above Normal")) %>%
  tab_source_note(
    source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Reference Climatolgy Period: 1982-2011.*"))
Average Temperature Anomalies - 2021
Degree Celsius Above Normal
Region Jan Feb Mar Apr May Jun
Global Oceans 0.24 0.24 0.28 0.30 0.31 0.32
Northeastern U.S. Shelf 1.20 0.91 0.91 0.93 1.47 1.74
Gulf of Maine 1.80 1.63 1.01 1.20 1.84 2.43
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Reference Climatolgy Period: 1982-2011.

Heatwaves Days

# 2. HW Days
heat_avgs %>% 
  gt(rowname_col = "Region") %>% 
  tab_stubhead(label = "Region") %>% 
  tab_header(
    title = md("**Number of Heatwave Days - 2021**")) %>%
  tab_source_note(md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Reference Climatolgy Period: 1982-2011.*"))
Number of Heatwave Days - 2021
Region Jan Feb Mar Apr May Jun
Global Oceans 31 28 31 30 31 28
Northeastern U.S. Shelf 17 14 15 15 20 22
Gulf of Maine 29 28 22 24 31 28
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Reference Climatolgy Period: 1982-2011.

Warming Rates Comparison

For each area the annual increase in temperature has been calculated for the period of 1982-2020 to account for all complete years of data available with OISSTv2.

Warming Rates

# 4. Get warming rates
area_rates <- map(area_hw, function(x){
  x <- x %>% 
    mutate(yr = year(time)) %>% 
    filter(between(yr, 1982, 2020)) %>% 
    group_by(yr) %>% 
    summarise(avg_temp = mean(sst, na.rm = T),
              avg_anom = mean(sst_anom, na.rm = T))
  
  temp_lm <- lm(avg_temp ~ yr, data = x)
  anom_lm <- lm(avg_anom ~ yr, data = x)
  return(
    list(
      temp = temp_lm,
      anom = anom_lm
    )
  )
  
})


# 5. Pull Rates
rate_table <- imap_dfr(area_rates, function(mod, area){
  tibble("Region"        = rep(area, 2),
         "Unit"          = c("Sea Surface Temperature", "Temperature Anomaly"),
         #"Intercept"     = c(coef(mod$temp)[[1]], coef(mod$anom)[[1]]),
         "Annual Change" = c(coef(mod$temp)[[2]], coef(mod$anom)[[2]])) %>% 
    mutate(`Annual Change` = round(`Annual Change`, 3))
})

Rate Table

rate_table %>% 
  filter(Unit == "Sea Surface Temperature") %>% 
  select(-Unit) %>% 
  gt(rowname_col = "Region") %>% 
  tab_stubhead(label = "Region") %>% 
  tab_header(
    title = md("**Annual Change in Sea Surface Temperature**"),
    subtitle = "Data from Years 1982-2020") %>%
  tab_source_note(md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Reference Climatolgy Period: 1982-2011.*"))
Annual Change in Sea Surface Temperature
Data from Years 1982-2020
Region Annual Change
Global Oceans 0.013
Northeastern U.S. Shelf 0.035
Gulf of Maine 0.041
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Reference Climatolgy Period: 1982-2011.

Heatwave Timelines

# Plotting Function
plot_mhw <- function(timeseries_data){
  
  
  # Set colors by name
  color_vals <- c(
    "Sea Surface Temperature" = "royalblue",
    "Heatwave Event"          = "darkred",
    "Cold Spell Event"        = "lightblue",
    "MHW Threshold"           = "coral3",
    #"MHW Threshold"           = "gray30",
    "MCS Threshold"           = "skyblue",
    #"MCS Threshold"           = "gray30",
    "Daily Climatology"        = "gray30")
  
  
  # Set the label with degree symbol
  ylab <- expression("Sea Surface Temperature"~degree~C)
  
  
  
  # Plot the last 365 days
  p1 <- ggplot(timeseries_data, aes(x = time)) +
      geom_segment(aes(x = time, xend = time, y = seas, yend = sst), 
                   color = "royalblue", alpha = 0.25) +
      geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), 
                   color = "darkred", alpha = 0.25) +
      geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
      geom_line(aes(y = hwe, color = "Heatwave Event")) +
      geom_line(aes(y = cse, color = "Cold Spell Event")) +
      geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = 1) +
      geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = 1) +
      geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = 1) +
      scale_color_manual(values = color_vals) +
      scale_x_date(date_labels = "%b", date_breaks = "1 month") +
      theme(legend.title = element_blank(),
            legend.position = "top") +
      labs(x = "", 
           y = ylab, 
           caption = paste0("Climate reference period :  1982-2011"))
    
  
  return(p1)
}

World

area_hw$`Global Oceans` %>% 
  filter(year(time) == 2021) %>% 
  plot_mhw()

Northeastern US Shelf

area_hw$`Northeastern U.S. Shelf` %>% 
  filter(year(time) == 2021) %>% 
  plot_mhw()

Gulf of Maine

area_hw$`Gulf of Maine` %>% 
  filter(year(time) == 2021) %>% 
  plot_mhw()

How Hot is 2021

This section seeks to place the monthly temperatures in the context of previous years for the purpose of answering:

  • Is this the hottest May on record?
  • If not, when was it hotter?
  • If so, when was the runner up?
  • Top 5?

Summary Prep

Regional Maps

For the Gulf of Maine and the Northeast Shelf a specific polygon has been used to both outline the area of interest and define what data to use for any of the area-specific calculations. These have been overlayed on the maps in blue to signify what areas are which.

Polygons for Mapping

# Polygons for mapping
new_england <- ne_states("united states of america") %>% st_as_sf(crs = 4326) 
canada      <- ne_states("canada") %>% st_as_sf(crs = 4326)
world       <- ne_countries() %>% st_as_sf(crs = 4326)
greenland   <- ne_states(country = "greenland") %>% st_as_sfc(crs = 4326)

# Load the crop areas
gom_shape <- get_timeseries_paths("gmri_sst_focal_areas")[["apershing_gulf_of_maine"]]["shape_path"]
gom_shape <- read_sf(gom_shape)
nelme_shape <- get_timeseries_paths("nelme_regions")$NELME$shape_path
nelme_shape <- read_sf(nelme_shape)

Load Satellite Data

# What data we want to load
data_window <- data.frame(lon = c(-720, 720),
                          lat = c(-90, 90),
                          time = as.Date(c("2021-01-01", "2021-06-30")))

# Load it up
oisst_daily <- oisst_window_load(oisst_path = oisst_path, data_window = data_window, anomalies = TRUE)
oisst_daily <- raster::stack(oisst_daily)

# fix the extent that is jank for no reason
extent(oisst_daily) <- extent(matrix(data = c(0, -90, 360, 90), nrow = 2))

Make it Monthly

# Make it monthly
make_monthly <- function(daily_ras){
  # Months to subset with
  month_key <- str_pad(c(1:12), 2, "left", 0) %>% setNames(month.abb)

  # names to match index to
  layer_index <- names(daily_ras)
  month_index <- str_sub(layer_index, 7, 8)
  
  # Pull distinct months
  months_present <- unique(month_index)
  month_key <- month_key[which(month_key %in% months_present)]
  
  # Pull the indices that match, take means
  map(month_key, function(x){
    
    # Pull days in month
    days_in_month <- which(month_index == x)
    
    # Take mean of those days
    month_avg <- mean(daily_ras[[days_in_month]])
  }) %>% 
    setNames(names(month_key))
  
  }


# Make it monthly
oisst_monthly <- make_monthly(oisst_daily)
monthly_stars <- map(oisst_monthly, ~ st_as_stars(rotate(.x)))

Universal Map Settings

# Set Labels for all the Plots
# plot_month <- "May"           # For testing single months
plot_months <- month.abb[1:6] # For Animating through each
plot_months <- setNames(plot_months, plot_months)


# Color limit
temp_limits <- c(-5, 5)
#temp_limits <- max(abs(values(oisst_monthly[[plot_month]])), na.rm = T) * c(-1,1) # Dynamic Limits

# Avg Anomaly for midpoint
color_midpoint <- 0
# color_midpoint <- mean(values(oisst_monthly[[plot_month]]), na.rm = T) # Dynamic midpoint

World

# Plot global map
plot_world <- function(plot_month){
  
  # text formatting for labels
  plot_lab <-  str_c("Monthly Temperature Anomalies - ", plot_month, " - 2021")
  guide_lab <- expression("Temperature Anomaly"~~degree~C)
  
  
  p1 <- ggplot() +
    geom_stars(data = monthly_stars[[plot_month]]) +
    geom_sf(data = world, fill = "gray90", size = .25) +
    scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                         limit = temp_limits, oob = scales::squish) +
    guides("fill" = guide_colorbar(title = guide_lab, 
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
    coord_sf(expand = F) +
    map_theme +
    labs(subtitle = plot_lab)
  
  return(p1)
}

# build list
world_maps <- map(plot_months, plot_world)

Jan

world_maps$Jan

Feb

world_maps$Feb

Mar

world_maps$Mar

Apr

world_maps$Apr

May

world_maps$May

Jun

world_maps$Jun

Northeastern US Shelf

# Map of Northeastern US Shelf

plot_shelf <- function(plot_month){
  
  # text formatting for labels
  plot_lab <-  str_c("Monthly Temperature Anomalies - ", plot_month, " - 2021")
  guide_lab <- expression("Temperature Anomaly"~~degree~C)
  
 

  p1 <- ggplot() +
    geom_stars(data = monthly_stars[[plot_month]]) + 
    geom_sf(data = new_england, fill = "gray90", size = .25) +
    geom_sf(data = canada, fill = "gray90", size = .25) +
    geom_sf(data = greenland, fill = "gray90", size = .25) +
    geom_sf(data = nelme_shape, fill = "transparent", size = .5, 
            color = gmri_cols("gmri blue")) +
    scale_x_continuous(breaks = seq(-180, 180, 5)) +
    scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                         limit = temp_limits, oob = scales::squish) +
    guides("fill" = guide_colorbar(title = guide_lab, 
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
    map_theme +
    coord_sf(xlim = c(-76.5, -64), ylim = c(35, 45.25))  +
    labs(subtitle = plot_lab)


  return(p1)
}


# build list
shelf_maps <- map(plot_months, plot_shelf)

Jan

shelf_maps$Jan

Feb

shelf_maps$Feb

Mar

shelf_maps$Mar

Apr

shelf_maps$Apr

May

shelf_maps$May

Jun

shelf_maps$Jun

Gulf of Maine

# Map of Gulf of Maine

plot_gom <- function(plot_month){

  # text formatting for labels
  plot_lab <-  str_c("Monthly Temperature Anomalies - ", plot_month, " - 2021")
  guide_lab <- expression("Temperature Anomaly"~~degree~C)
  
  
  p1 <- ggplot() +
    geom_stars(data = monthly_stars[[plot_month]]) + 
    geom_sf(data = new_england, fill = "gray90", size = .25) +
    geom_sf(data = canada, fill = "gray90", size = .25) +
    geom_sf(data = greenland, fill = "gray90", size = .25) +
    geom_sf(data = gom_shape, fill = "transparent", size = .5, 
            color = gmri_cols("gmri blue")) +
    scale_x_continuous(breaks = seq(-180, 180, 5)) +
    scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                         limit = temp_limits, oob = scales::squish) +
    guides("fill" = guide_colorbar(title = guide_lab, 
                                   title.position = "top", 
                                   title.hjust = 0.5,
                                   barwidth = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
    map_theme +
    coord_sf(xlim = c(-72.5, -64), ylim = c(39, 45.25))  +
    labs(subtitle = plot_lab)

  return(p1)

}
# build list
gom_maps <- map(plot_months, plot_gom)

Jan

gom_maps$Jan

Feb

gom_maps$Feb

Mar

gom_maps$Mar

Apr

gom_maps$Apr

May

gom_maps$May

Jun

gom_maps$Jun